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FOREWORD 


This  manual  is  designed  to  be  used  by  participants  in  a  2-day  workshop  on 
the  MOPAC  computational  chemistry  program.  MOPAC  implements  the 
semi -empirical  molecular  orbital  methods  pioneered  by  Professor  Michael  J.  S. 
Dewar,  University  of  Texas  at  Austin.  This  manual  is  aimed  at  beginning 
users,  assuming  some  fundamental  knowledge  of  quantum  mechanics,  but  not 
requiring  familiarity  with  the  higher  mathematics  involved  in  molecular 
orbital  calculations.  This  second  edition  adds  three  appendices  containing 
information  which  was  found  to  be  of  interest  to  attendees  of  the  first  three 
workshops.  Only  the  material  in  Appendix  A  is  definitely  covered  in  the 
workshop.  Appendices  B  and  C  do  not  need  to  be  covered  to  achieve  the  goals 
of  the  workshop.  The  workshop  addresses  one  objective  in  a  joint  U.S.  Air 
Force  Academy  Department  of  Chemistry  —  Frank  J.  Seiler  Research  Laboratory 
project  to  establish  a  center-of -excellence  in  computational  chemistry  at  the 
Air  Force  Academy.  That  objective  is  to  promote  use  of  MOPAC  in  DOD  and 
associated  research  groups  while  demonstrating  its  applicability  to  real 
problems  of  interest  to  the  Air  Force.  The  project  is  funded  by  the  Air  Force 
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MOPAC  WORKSHOP 


Outline  of  workshop: 


1.  Overview  of  the  MOPAC  method. 


2.  Description  of  theoretical/computational  process. 


3.  How  to  use  MOPAC. 


4.  How  to  use  DRAW. 


5.  Examples  of  applications. 


6.  Other  considerations. 


1.  Overview  of  the  MOPAC  method. 

MOPAC  is  a  semi-empirical  molecular  orbital  program  for  calculating 
various  properties  of  molecular  and  ionic  systems.  The  geometries,  charge 
distributions,  dipoles,  the  make-up  of  the  molecular  orbitals,  and  other 
similar  information  on  systems  of  Interest  can  be  calculated.  Thermodynamic 
information  can  be  obtained  from  the  calculated  heats  of  formations  (AHf)  of 
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systems.  Kinetic  information  can  be  obtained  from  calculations  which  map  the 
change  in  AHf  with  chosen  reaction  coordinates  and  identify  the  transition 
states  for  reactions.  The  Infrared  spectra  of  molecules  can  be  calculated  for 
comparison  with  experimental  spectra.  Other  more  specialized  capabilities, 
such  as  the  calculation  of  polymer  properties,  are  also  available.  The  MOPAC 
Manual  is  the  document  which  shows  how  the  full  range  of  MOPAC' s  capabilities 
can  be  used.  This  workshop  is  intended  to  get  you  started  as  a  MOPAC  user  and 
to  give  you  a  sense  of  its  powers  and  limitations. 

As  the  term,  "semi-empirical ,"  implies,  a  MOPAC  calculation  is  based 
partly  on  quantum  mechanical  theory  and  partly  on  experimental  data. 

Theoretical  considerations  led  to  the  form  of  the  Hamiltonian  energy  operator 
used  In  the  secular  equations,  i he  secular  equations  are  solved  using  matrix 
algebra  to  obtain  the  system  eigenfunctions  ard  the  corresponding 
eigenvalues.  These  eigenvalues  are  then  used  to  calculate  the  system's 
electronic  energy  (EE).  EE  is  combined  with  the  calculated  nuclear  repulsion 
(EN),  the  energies  of  the  Isolated  atoms  (EISOL),  and  the  heats  of  atomization 
of  the  atoms  in  the  system  (EHEAT)  to  give  the  heat  of  formation  of  the 
system.  The  diagram  below  illustrates  this  process  for  the  example  of  acetone. 


AHf 


3C  (gr)  +  3H?  +  1/2  O2 


v  EHEAT 


3C  (g)  +  6H  *  0 


1 

24  e- 


EISOL 
+  3C4+ 


+ 


»  CH3COCH3 
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There  are  three  different  Hamiltonians  which  were  developed  by  the  research 
group  of  Michael  J.  S.  Dewar  at  the  University  of  Texas  at  Austin  and  are 
available  In  the  MOPAC  program:  MNDO,  AMI,  and  MINDO/3. 

Experimental  data  enter  into  the  method  because  the  secular  equations  set 
up  from  these  Hamiltonians  contain  functions  (describing  the  behavior  of  each 
atom  In  the  system)  which  Include  adjustable  parameters.  These  parameters  are 
determined  not  from  first  principles  (i.e.,  purely  theoretical  considera¬ 
tions),  but  rather  by  a  fitting  process  called  "parameterization."  Parameter¬ 
ization  uses  experimental  data  from  a  set  of  molecules  to  determine  the  values 
of  these  parameters.  There  are  seven  parameters  for  each  element  in  the  MNDO 
method.  The  AMI  method  has  the  same  seven  and  a  few  more  per  element  while 
the  MINDO/3  method  has  seven  per  pair  of  atoms.  In  the  fourth  edition  of 
MOPAC,  the  atoms  which  have  been  parameterized  are:  H,  Li,  B,  C,  N,  0,  F,  Al, 
SI,  P,  S,  Cl,  Ge,  Br,  Sn,  I,  Hg,  Pb.  (Currently  only  MNDO  is  parameterized 
for  all  these  atoms.) 

MOPAC  can  be  run  on  normal  Control  Data  Corporation  (CDC),  Data  General, 
Gould,  and  Digital  Equipment  Corporation  (VAX)  computers,  as  well  as  on  the 
CDC  205  and  Cray-XMP  "supercomputers."  Versions  are  also  available  for 
several  microcomputers , including  IBM  PC-AT  and  XT  and  Zenith. 

The  process  In  which  MOPAC  Is  a  valuable  tool  begins  with  a  clear 
statement  of  the  chemical  problem  being  dealt  with:  Is  a  particular  molecule 
stable?  What  is  the  most  stable  geometry  of  a  particular  ion?  What  is  the 
electronic  structure  of  a  molecule?  What  is  the  enthalpy  of  reartion  for  a 
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particular  reaction?  How  large  is  the  energy  of  activation  for  the  reaction 
when  it  proceeds  in  a  certain  way?  What  are  the  frequencies  of  the  normal 
modes  of  a  particular  molecule?  Which  compound  requires  the  least  amount  of 
activation  to  react  with  a  specific  substrate?  And  so  on. 

Chemists  answer  these  kinds  of  questions  every  day  without  MOPAC.  They 
get  their  answers  from  thermodynamic  tables,  the  literature,  or  experiments 
they  perform  themselves.  MOPAC  is  useful  when  these  sources  cannot  easily  or 
reliably  provide  the  desired  information.  But  its  greatest  value  emerges  when 
it  is  used  to  predict,  interpret,  and  confirm  experimental  results. 

Once  the  system  and  the  Information  desired  about  it  are  determined,  the 
MOPAC  user  constructs  a  data  file  which  models  the  :ystem  and  describes  the 
desired  information  in  a  form  which  the  MOPAC  program  can  use  to  carry  out  the 
needed  calculations.  The  essential  parts  of  the  data  file  are  the  KEYWORDS, 
which  tell  the  program  how  the  calculation  is  to  be  carried  out  and  what 
output  Is  desired,  and  the  GEOMETRY  DEFINITION,  which  gives  the  coordinates  of 
the  atoms  in  the  system.  The  user  then  enters  a  simple  command  to  start  the 
calculation  and  the  MOPAC  program  takes  over.  It  performs  all  the  requested 
calculations  and  creates  a  number  of  files  containing  the  results.  Most  of 
the  desired  Information  will  be  in  an  output  file.  In  the  third  section  of 
this  manual,  we  will  show  how  to  create  a  data  file  (labeled  FILENAME.DAT)  and 
how  to  interpret  an  output  file  (labeled  FILENAME. OUT) .  But  first  we  will 
give  in  the  following  section  a  brief  outline  of  what  goes  on  when,  as  we  said 
above,  "the  MOPAC  program  takes  over." 
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2.  Description  of  theoretical/computational  process. 

First  MOPAC  check,  the  KEYWORDS  for  compatibility.  For  example.  If  both 
AMI  and  MINDO-3  were  entered,  the  calculation  would  stop  since  only  one 
Hamiltonian  can  be  used.  It  also  checks  that  the  geometry  definition  is  not 
fatally  flawed,  e.g.,  two  atoms  positioned  at  the  same  point  in  space.  If  the 
data  file  passes  this  inspection,  the  subroutine  MOIDAT  is  called. 

MOLDAT  retrieves  (from  a  set  of  reference  data)  the  parameters  and  other 
atomic  data  for  the  atoms  comprising  the  system.  It  sets  up  the  many  arrays 
(or  multi-dimensional  matrices)  to  be  used  In  the  calculations  based  on  the 
total  number  of  electrons  In  the  system  and  other  factors,  such  as  whether  a 
UHF  (unrestricted  Hartree-Fock)  or  an  RHF  (restricted  Hartree-Fock ,  this  being 
the  default  selection)  calculation  has  been  requested. 

MOPAC  then  proceeds  to  one  of  the  seven  subroutines  shown  in  the  flow 
diagram  in  Fig  1  depending  on  which  KEYWORD  has  been  selected.  We  will  follow 
what  happens  when  none  of  them  have  been  selected  so  that  the  default 
subroutine,  FLEPO,  Is  called.  (This  is  a  geometry  optimization  procedure 
still  named  after  two  of  its  creators,  Fletcher  and  Powell,  although  with  the 
fourth  edition,  another  procedure,  the  BFGS  method,  is  used).  The  GEOMETRY 
package  (of  which  FLEPO  Is  a  part)  communicates  with  the  .ELECTRONICS  package 
through  subroutine,  COMPFG.  These  two  packages  comprise  about  90%  of  the 
program.  Following  the  basic  workings  of  FLEPO  requires  some  understanding  of 
the  terms  listed  below.  In  this  list  and  in  the  following  discussion,  a 
singly  underlined  symbol  indicates  a  vector  and  a  doubly  underlined  symbol 
indicates  a  matrix. 
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Fig.  1.  MOPAC  flow  diagram. 


xj°  -  the  1th  coordinate  of  the  3N-6  coordinates  (where  N  is  the  number  of 
atoms  in  the  system)  initially  entered  in  the  geometry  definition  of  the  data 
file. 


Xjk  -  the  1th  coordinate  after  completion  of  the  kth  cycle  of  the  calculation. 
What  a  cycle  consists  of  will  become  clear  as  we  proceed. 

xk  -  the  vector  in  <3N-6)-space  with  the  internal  coordinates  as  components 
(x1k,x2k,X3k,...,x(3N_6)k) 

g^k  -  the  derivative  of  AHf  with  respect  to  in  the  kth  cycle  of  the 
calculation. 

-  the  vector  with  the  gjk  as  components,  i.e.,  (gjk,  g2k,  93k. • •  g(3N-6)k>- 
pk  -  the  change  in  x  from  one  cycle  to  the  next,  i.e.,  pk  =  xk+1  -  xk. 

yk  -  the  change  in  3  from  one  cycle  to  the  next,  i.e.,  yk  =  c[k+1  - 

-  the  Hessian  matrix,  this  is  the  matrix  of  the  second  derivatives  of  AHf 
with  respect  to  the  coordinates,  d2AHf /dx^dxj ,  on  the  kth  cycle.  Thus  it  is  a 
square  matrix  with  3N-6  rows. 

-  the  inverse  Hessian,  i.e.,  the  Inverse  of  £k. 

dk  -  the  "search  direction  vector"  used  to  calculate  gk.  It  is  given  by  dk  = 

akak 
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a  -  a  scaling  factor  for  which  accounts  for  the  fact  that  the  approach  to 
minimum  energy  on  the  3N-6  dimensional  surface  will  not  be  along  a  pure 
parabolic  surface. 

With  these  definitions  In  hand,  let  us  follow  what  FLEPO  does.  Using  the 
coordinates,  x,0,  which  you  entered  in  the  data  file,  it  performs  a 
self-consistent  field  molecular  orbital  (SCF  MO)  calculation  in  the 
ELECTRONICS  package  using  the  MNOO,  AMI,  or  MINDO/3  methods.  We  will  not  go 
into  the  details  of  how  these  calculations  are  done  since  they  are  described 
in  their  essential  form  in  any  introductory  quantum  mechanics  text.*  The 
result  is  the  initial  AHf°  for  the  system.  FLEPO  also  calculates  the  values 
for  gj°  in  the  GEOMETRY  package.  The  coordinates,  x^,  to  be  used  in  the  next 
cycle  are  then  calculated  from 

x \ 1  =  xi°  -  0.01*  SIGN(gi °> , 

where  SIGN(g^°)  »  1  if  g^0  >  0  and  =  -1  otherwise.  The  second  term  in  this 
equation  corresponds  to  the  1th  component  of  p*.  which,  In  subsequent 
calculations,  will  be  calculated  differently,  as  we  shall  see.  With  the  new 
values,  x^,  the  next  cycle,  which  also  consists  of  only  a  single  SCF 
calculation  (subsequent  cycles  typically  consist  of  four  SCF's),  is  completed 
and  the  results,  AHf1  and  g1 ,  are  obtained.  The  vectors,  (y°  =  g1  -  g°)  and 


*also  see  "MOPAC:  A  Semi-Empirical  Molecular  Orbital  Program,"  James  J.  P. 
Stewart,  Handbook  of  Computer  Aided  Modeling  and  Design,  ACS  Monograph  Series, 
1988. 
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(E°  *  x1  -  x°>.  are  then  used  to  calculate  the  diagonal  elements,  Gi,  = 
yi°/pi°,  of  the  Hessian  matrix,  so  that  ^  and  d1*-  can  be  obtained  (for 
k=l)  as  defined  above.  FLEPO  next  calls  the  subroutine,  LINMIN,  to  calculate 
a  value  for  a  using  a  parabolic  fit.  With  a  and  d^  in  hand,  LINMIN  calculates 
from 


=  adk  . 

Now  FLEPO  can  calculate  x  for  the  next  cycle  (k+l=2)  using 

xk+1  = 

The  next  cycle  results  in  new  values  for  AHf ,  g,  £,  d,  and  a,  as  well  as  the 
coordinates  for  the  following  cycle  using  the  same  equations.  The  equation 
used  to  calculate  the  inverse  of  the  Hessian  in  a  following  cycle  is  given 
below: 


yk+1  _  jjk  -(flkykgkt  +  gkyktjjk)/5  +  Q(gkgkt>/s  t 

where  0=1+  y^tjjkyk/gktyk 
and  S  = 

(Note:  the  superscript,  t,  Indicates  the  transpose  of  the  vector) 

Hhat  must  be  kept  in  mind  at  this  point  is  that  the  Inverse  Hessian  and 
all  the  other  quantities  described  above  are  calculated  to  obtain 
progressively  better  coordinates,  xj,  for  the  system.  Coordinates  get 
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"better"  when  the  heat  of  formation  calculated  from  them  is  lower  than  that 
obtained  in  the  previous  cycle. 

After  each  cycle  FLEPO  calculates  an  indicator  of  how  close  the  calculated 
value  of  AHf  is  to  a  minimum.  This  indicator,  GNORM,  is  given  by 

GNORM  *  / 1  g^2  =  |j[|  . 

which  has  the  units  kcal/mol/A  and  can  be  considered  a  measure  of  how 
sensitive  AHf  is  to  small  changes  In  the  system  geometry. 

FLEPO  continues  to  run  through  these  cycles  until  two  tests  are  passed: 

(1)  the  |gi|  values  (i.e.,  the  individua1  components  of  GNORM)  must  be 
below  a  certain  limit  which  is  a  calculated  function  of  the  system,  and 

(2)  any  one  of  the  following: 

(a)  "test-on-x-sati sf led" ,  i.e.,  the  magnitude  of  the  vector,  d*.  is 
very  small  (less  than  0.0001  A), 

(b)  Herbert's  test,  i.e.,  the  projected  decrease  in  AHf  is  less  than 
0.001  kcal /mol , 

(c)  "gradient  satisfied",  i.e.,  GNORM  is  less  than  KEYWORD-specif ied 
value  (default  is  1  kcal/mol/A), 
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(d)  "AHf  test  satisfied",  i.e.,  the  change  in  AHf  from  cycle  to 
cycle  is  less  than  0.002  kcal/mol. 

While  the  process  we  have  described  is  controlled  by  only  one  of  the 
subroutines  shown  in  the  flow  diagram,  it  is  a  very  important  process  which 
occurs  in  most  of  the  MOPAC  calculations  which  are  done.  Brief  descriptions 
of  what  is  done  when  KEYWORDS  select  the  other  subroutines  are  given  in  the 
MOPAC  Manual.  The  description  above  was  given  not  only  to  satisfy  the 
curiosity  of  users  who  prefer  to  have  some  understanding  of  what  happens  when 
"the  MOPAC  progam  takes  over,"  but  also  to  introduce  some  terms  which  must  be 
understood  to  properly  interpret  MOPAC  results. 

3.  How  to  use  MOPAC. 

We'll  start  by  calculating  the  heat  of  formation  of  a  simple,  well  known 
molecule,  methanol  (so  a  good  filename  would  be  CH30H).  We  need  a  computer  on 
which  MOPAC  has  been  loaded  and  access  to  a  directory  in  which  we  can  create  a 
data  file  and  initiate  the  MOPAC  calculation.  We  then  enter  a  command  to  edit 
the  data  file,  e.g.,  EDIT/EDT  CH30H.DAT  for  a  VAX  or  vi  ch3oh.dat  for  a  UNIX 
operating  system.  Since  that  file  presumably  does  not  exist,  we  get  a  blank 
screen  and  begin  creating! 

a.  Creating  the  data  file,  CH30H.DAT. 

The  first  line  will  be  a  series  of  KEYWORDS,  which  tell  the  program  which 
properties  we  want  calculated  for  the  system  we  will  enter,  how  much  CPU  time 
it  will  be  allowed,  the  precision  we  want,  what  outputs  we  want,  etc.  The 
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full  list  of  KEYWORDS  and  their  effects  on  the  calculation  is  given  in  Chapter 
2  of  the  Manual.  In  every  calculation  a  AHf,  geometry,  and  charge  di st i but i on 
for  the  system  are  output  whether  you  want  them  or  not.  So  for  our  simple 
calculation,  this  first  line  will  look  like  this: 

T=3600  AMI  PRECISE 

T=3600  gives  the  computer  3600  seconds  of  CPU  time  to  complete  the 
calculation,  AMI  tells  the  program  to  use  the  AMI  method  rather  than  the  MNDO 
method,  which  is  the  default  choice,  and  PRECISE  puts  into  effect  criteria 
which  will  give  us  a  more  precise  result  than  the  default  criteria,  but  at  the 
expense  of  more  computer  time. 

The  second  and  third  lines  are  comments  which  we  may  use  to  describe  for 
ourselves  what  the  calculation  is,  for  example: 

AMI  CALCULATION  OF  HEAT  OF  FORMATION  OF  METHANOL  WITH  PRECISE 

The  first  of  these  three  lines  must  always  be  used  for  the  purpose 
described  above,  and  the  second  and  third  lines  must  be  left  blank  if  no 
comments  are  desired. 

Now  we  get  to  the  starting  geometry  definition.  The  geometry  is  defined 
in  internal  coordinates  with  a  connectivity  list,  which  is  the  list  of  numbers 
in  three  columns  under  the  heading  "Connectivity"  in  the  table  below.  How  it 
works  should  become  clear  as  we  discuss  this  table.  The  numbers  in  the 
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columns,  "Distance,"  "Angle,"  and  "Dihedral"  are  the  3N-6  coordinates,  xj°, 
which  were  discussed  in  Section  2.  The  same  geometry  can  be  defined  in  many 
ways,  and  what  is  shown  below  is  just  one  example: 


Atom 

Distance 

Angle 

Dihedral 

Connectivity 

I 

I-A 

I-A-B 

I-A-B-C 

A 

B  C 

1 

i 

iC 

2 

to 

1.43  1 

0  0 

0 

0 

1 

3 

IH 

0.96  1 

108.0  1 

0 

0 

2 

1 

4 

lH 

1.10  1 

109.0  1 

180.0 

1 

1 

2  3 

5 

H 

1.10  1 

109.0  1 

-60.0 

1 

1 

2  3 

6 

'h 

1 

1.10  1 

109.0  1 

BLANK  LINE.. 

60.0 

1 

1 

2  3 

The  system  described  above  is  shown  in  Fig  2.  The  drawing  was  created  from 
the  CH30H.DAT  file  using  the  DRAW  program,  which  will  be  discussed  later.  The 
geometry  definition  can  be  understood  by  following  the  definition  for  one 
atom,  say  atom  number  6,  a  hydrogen  designated  H(6). 

Distance  I-A:  H(6)  is  1.10  A  from  C( 1 )  since  in  the  connectivity  list,  A  is 
defined  to  be  atom  number  1.  This  is  the  distance  6-1. 

Angle  I-A-B:  H(6)  forms  an  angle  of  109.0  degrees  with  C(l)  and  0(2)  since  in 
the  connectivity  list,  A  and  B  are  defined  to  be  atoms  1  and  2,  respectively. 
This  is  the  angle  6-1-2. 

Dihedral  I-A-B-C:  The  H(6)-C(l)  bond  is  rotated  60.0  degrees  clockwise  from 
the  0(2)-H(3)  bond  when  one  looks  down  the  0(2)-C(l)  bond,  with  the  0(2) 
closer  to  the  eye.  See  Fig  2.  This  is  the  dihedral  6-1-2-3. 
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Fig.  2.  CH30H  geometry. 
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The  number  "1"  following  a  coordinate  tells  MOPAC  to  optimize  that 
coordinate.  A  "0"  in  this  position  tells  MOPAC  to  keep  that  coordinate 
fixed.  There  must  be  at  least  one  space  between  each  entry  in  a  row,  but  we 
normally  put  in  more  for  ease  in  reading  the  data  file.  The  0's  in  the  rows 
for  0(2)  and  K'3)  are  used  to  indicate  that  the  entries  can  not  be  defined  and 
are  there  only  to  maintain  format/spacing  requirements  up  to  the  connectivity 
list. 

The  atom  numbers  and  the  column  headings  shown  in  the  geometry  description 
above  cannot  be  put  in  the  data  file.  Only  the  part  within  the  dashed  line  is 
entered.  Thus  the  completed  CH30H.DAT  file  will  look  like  this: 

T-3600  AMI  PRECISE 

AMI  CALCULATION  OF  HEAT  OF  FORMATION  OF  METHANOL  WITH  PRECISE 
C 


0 

1  .43 

1 

0 

0 

0 

0 

1 

H 

0.96 

1 

108.0 

1 

0 

0 

2 

1 

H 

1.10 

1 

109.0 

1 

180.0 

1 

1 

2 

3 

H 

1.10 

1 

109.0 

1 

-60.0 

1 

1 

2 

3 

H 

1.10 

1 

109.0 

1 

60.0 

1 

1 

2 

3 

(NOT  FORGETTING  THAT  THE  KEYWORD  LINE  IS  THE  FIRST  LINE  IN  THE  FILE)  Now  a 
command  is  entered  to  start  the  calculation  or  at  least  submit  it  to  a  queue 
where  it  will  wait  its  turn  (e.g.,  MOPAC  CH30H  60  on  a  VAX,  where  the  60  is  an 
indicator  of  the  length  of  time  estimated  for  the  job  and  thus  the  queue  to 
which  the  job  should  be  sent). 

Next  you  go  for  a  cup  of  coffee,  lunch,  or  a  holiday,  depending  on  how 
many  atoms  are  in  the  system,  the  kind  of  computer  you're  using,  and  how  busy 


15 


it  is.  With  CH30H  using  a  not-too-busy  VAX,  you'll  probably  have  just  enough 
time  to  pour  yourself  a  cup  of  coffee  before  the  job  is  done.  Assuming  it  has 
been  successfully  completed,  you  should  now  have  some  additional  files  in  your 
directory,  e.g.,  CH30H.0UT  and  CH30H.ARC  on  the  VAX.  Next  let's  take  a  look 
at  how  the  results  of  the  calculation  are  presented  and  what  they  mean. 

b.  Reading  the  output  file,  CH30H.0UT. 

It  must  be  kept  in  mind  that  while  the  kind  of  output  and  the  order  in 
which  it  is  presented  in  this  file  is  fairly  typical,  changes  in  keywords  can 
make  an  output  file  look  quite  different.  Understanding  the  parts  of  the  file 
we  will  discuss  below  will  be  very  helpful  in  reading  any  MOPAC  output  file. 

(1)  After  some  initial  labeling,  the  KEYWORDS  are  listed.  These 
should  be  checked  to  Insure  that  this  is  the  result  of  the  calculation  you 
think  it  is! 

(2)  Next  the  essential  ingredients  of  the  data  file  used  are 
presented.  This  listing,  especially  the  geometry  definition,  lets  you  know  on 
what  system  this  calculation  was  done. 

(3)  Following  the  listing  of  the  CARTESIAN  COORDINATES  is  a 
statement  of  how  many  singly  and  doubly  occupied  energy  levels  (or  MO's)  were 
calculated  for  the  system.  This  can  be  revealing  if  you  expected  the  system 
to  be  a  singlet  in  the  ground  state  and  find  that  there  is  a  singly  occupied 
MO. 
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(4)  Next  come  the  interatomic  distances  for  your  initial  geometry 


followed  by  some  information  on  how  many  cycles  were  required  before  the  tests 
were  passed  which  signalled  the  calculation  to  end.  The  GNORM  (GRAD)  and  AHf 
(HEAT)  at  the  end  of  each  cycle  are  also  given. 

(5)  Then  come  the  main  results  of  the  calculation,  the  FINAL  HEAT  OF 
FORMATION  and  IONIZATION  POTENTIAL,  both  of  which  can  be  readily  compared  with 
experimental  results.  There  is  also  some  other  self-explanatory  Information 
given  on  the  system  and  the  calculation.  The  calculated  value  of  -57.0 
kcal/mol  does  not  agree  very  well  with  the  experimental  value  of  -48.2.  The 
error  of  -8.8  kcal/mol  is  larger  than  the  average  absolute  error  in  AHf 
obtained  using  the  methods  Included  in  MOPAC.  A  brief  table  giving  some 
examples  of  these  errors  Is  given  below. 

Examples  of  Average  Errors 

_ AM] _  _ MNDO _ 

Species  AHf  bond  bond  AHf  bond  bond 


Containing: 

(kcal /mol ) 

lenqth.  A 

(kcal /mol ) 

lenqth,  A  anqle 

CHNO* 

5.5 

0.014 

3° 

6.3 

0.014  2.8° 

CH 

5.1 

6.0 

F 

9.8 

Cl 

5.6 

Br 

7.2 

I 

5.8 

S 

7.7 

*138  compound  test  set. 
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(6)  Next  is  the  listing  of  the  internal  coordinates  for  the  final 
optimized  geometry.  Rather  than  studying  this  listing,  however,  we  normally 
use  the  DRAW  program  discussed  in  Section  4  to  examine  the  geometry.  Again  it 
is  Important  to  know  how  well  MOPAC  methods  predict  geometries.  The  table 
above  indicates  the  kinds  of  errors  in  geometries  you  can  expect  using  MOPAC 
methods . 

(7)  The  next  listing,  INTERATOMIC  DISTANCES,  gives  calculated 
distances  between  the  atoms  in  the  final  optimized  geometry.  Here  again  DRAW 
is  usually  a  more  desirable  vehicle  for  studying  this  aspect  of  the  system 
because  of  the  visual  displays  it  presents. 

(8)  The  energies  (in  eV>  of  the  MO  energy  levels  are  then  printed 
under  the  heading,  EIGENVALUES.  Each  value  corresponds  to  an  MO  which 
contains  at  most  two  electrons.  Thus  for  degenerate  orbitals,  two  values 
which  are  virtually  identical  will  appear  in  succession.  Since  the  number  of 
filled  levels  is  given  (See  (3)  above),  the  gap  between  the  highest  occupied 
MO  (HOMO)  and  lowest  unoccupied  MO  (LUMO)  can  be  determined  and  related  to 
reactivity  or  other  properties  of  the  system.  Note  that  these  are  the  MO's 
formed  only  from  the  valence  electrons  of  the  atoms  in  the  system. 

(9)  The  ATOM  ELECTRON  OENSITY  given  in  the  next  table  is  obtained 
for  each  atom  as  follows.  The  table  of  ATOMIC  ORBITAL  ELECTRON  POPULATIONS 
(two  tables  further  down  in  the  CH30H.0UT  file),  gives  the  electron  population 
of  each  atomic  orbital  in  the  basis  set  used  to  construct  the  LCAO  MO's.  They 
are  given  in  the  same  order  in  which  the  atoms  in  the  system  are  numbered. 
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Thus  the  first  4  values  give  the  electron  population  for  the  2s,  2px,  2py,  and 
2pz  orbitals  of  C(1).  (Only  the  valence  atomic  orbitals  are  Included.)  These 
values  are  the  diagonal  elements  of  the  DENSITY  MATRIX,  which  Is  essentially 
the  square  of  the  MO's  obtained  from  the  SCF  calculation.  The  sum  of  these 
four  electron  populations  gives  the  ATOM  ELECTRON  DENSITY  for  C(l),  4.0733  in 
the  CH30H  file.  Since  these  four  orbitals  would  carry  a  charge  of  4.0000  in 
the  unbonded  carbon,  the  CHARGE  column  shows  a  net  charge  on  C( 1 )  of  -0.0733. 
The  same  analysis  for  the  other  atoms  In  the  molecule  gives  the  overall  charge 
distribution  of  the  system. 

(10)  While  the  charge  distribution  cannot  be  directly  compared  with 
experimental  data,  the  total  dipole  moment  of  the  molecule  can  be  calculated 
from  these  results  and  compared  with  the  experimental  result.  It  is  given  in 
the  lower  right  corner  of  the  table  below  the  one  on  NET  ATOMIC  CHARGES  AND 
DIPOLE  CONTRIBUTIONS.  For  this  system  it  is  1.621  debyes. 

(11)  Finally  there  is  the  table  showing  the  optimized  geometry  in 
cartesian  coordinates. 

The  CH30H.ARC  file  is  an  abbreviated  version  of  the  .OUT  file  which  is 
convenient  for  long  term  storage  of  results.  It  is  also  the  file  which  is 
normally  used  by  the  DRAW  program  to  display  the  results  of  a  MOPAC 
calculation  in  graphic  form.  We  now  discuss  how  DRAW  can  be  used. 

4.  How  to  use  DRAW. 

DRAW  is  a  program  designed  to  give  graphical  displays  of  systems  whose 
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geometry  definitions  are  contained  in  output  and  data  files  in  MOPAC.  The 
DRAW  Manual  gives  a  complete  description  of  the  many  capabilities  of  this 
program.  Within  the  software  itself  there  is  a  HELP  routine  which  can  be 
called  to  assist  you  while  you  are  using  it.  In  this  workshop  we  will 
demonstrate  some  of  the  more  commonly  used  functions  of  DRAW. 

a.  Drawing  the  initial  geometry. 

Drawing  the  initial  geometry  in  a  data  file  is  almost  always  done  to 
insure  that  the  GEOMETRY  DEFINITION  is  what  you  intended.  After  exit’ng  the 
data  file,  you  issue  the  command.  DRAW,  to  get  into  the  program.  In  response 
to  a  prompt,  you  enter  the  filename,  CH30H.DAT,  to  bring  that  file  into  the 
program.  Next,  the  letter  "T"  followed  by  the  name  for  the  graphics  terminal 
you  are  using  is  entered,  e.g.,  T  4107  for  a  Tektronix  4107  terminal.  This 
produces  a  picture  of  the  system  in  the  data  file  using  a  stick  model  with 
atoms  represented  by  the  numbers  they  were  assigned  in  the  file.  Now  there 
are  a  number  of  ways  you  can  manipulate  this  oicture  to  get  a  better  look  at 
the  geometry.  Some  commonly  used  examples  are  given  below. 

(1)  You  can  rotate  the  molecule  by  issuing  a  command  such  as  ROTATE 
CARTESIAN  2  0  90  45,  which  will  rotate  the  molecule  around  the  number  2  atom 
in  cartesian  coordinates  0,  90,  and  45  degrees  about  the  x,  y,  and  z  axes, 
respectively.  Or  you  can  issue  the  command,  ROTATE  PAIR  1  2  90,  which  will 
rotate  the  system  90  degrees  about  the  axis  connecting  the  pair  of  atoms,  1-2. 

(2)  You  can  enter  the  DISPLAY  mode  (entering  DISPLAY)  and  change  the 
style  of  the  picture.  For  example,  a  popular  ball  and  stick  picture 
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designated  ORTEP  can  be  produced  by  entering  DISPLAY  ORTEP  QUIT,  where  QUIT 
tells  the  program  you  want  it  to  quit  the  display  mode  and  execute  the  change 
(to  ORTEP).  The  HELP  routine  gives  other  styles  which  are  available. 

(3)  You  can  change  the  way  the  atoms  are  labelled  by  entering 
DISPLAY  LABEL,  which  results  in  a  prompt  asking  you  whether  you  want  atom 
labels  to  be  their  chemical  symbols  (enter  SYMBOL),  removed  altogether  (enter 
MASK),  to  be  user-defined  (enter  USER),  etc.  Again,  once  you  have  specified 
your  choice,  enter  QUIT  to  put  the  change  into  effect. 

(4)  Bonds  are  drawn  between  atoms  which  are  close  enough  to  be 

within  the  sum  of  their  Van  der  Waals  radii  from  each  other.  This  criterion 
can  be  altered  by  entering  DISPLAY  BONDS,  which  results  in  a  prompt  for  you  to 
select  a  choice  of  different  criteria,  including  simply  removing  bonds  which 

you  do  not  wish  to  appear.  Again  QUIT  puts  the  changes  into  effect. 

(5)  With  the  system  drawn  on  the  terminal  you  can  get  information 

about  the  geometry  by  entering  INFORMATION  to  get  into  the  INFORMATION  mode. 

Then,  for  example,  you  can  find  the  distance  between  any  two  atoms.  Entering 
2  6  immediately  gives  the  distance  between  atoms  2  and  6.  Entering  236 
gives  the  angle  between  2,  3,  and  6.  Entering  1236  gives  the  dihedral 
angle. 

b.  Drawing  the  optimized  geometry. 

Drawing  the  optimized  geometry  following  a  successful  calculation  is 
obviously  the  best  way  to  see  how  it  turned  out.  You  generate  the  picture  in 
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the  way  described  above  for  the  data  file,  here  using  the  .ARC  or  .OUT  files. 
It  can  then  be  manipulated  in  the  same  way  as  described  above  for  the  data 
file  picture.  But  there  are  other  things  you  can  do  with  an  .OUT  or  .ARC 
file,  some  examples  of  which  are  given  below. 

(1)  Entering  the  INFORMATION  mode  with  the  command  I  I  0  gives  a 

summary  of  the  information  contained  in  an  ARC  file.  Or  you  can  ask.  for  other 

specific  information  about  the  system  as  outlined  in  the  HELP  routine. 

(2)  If  the  ARC  file  contains  several  optimized  geometries  because  a 
reaction  path  calculation  was  done  (more  on  this  in  the  next  section),  then 
the  various  points  along  the  reaction  path  can  be  viewed  in  sequence  by 
entering  NEXT  after  each  one  is  shown. 

(3)  If  the  KEYWORD,  "FORCE,"  was  used  to  get  the  vibrational  modes 

of  the  molecule  (more  on  this  later),  then  the  .OUT  file  must  be  brought  into 

DRAW.  Now  the  various  modes  can  be  viewed  by  entering  a  command  such  as 

DISPLAY  VIBRATION  14  QUIT,  which  results  in  the  molecule  being  displayed  in 
stick  style  with  arrows  indicating  the  direction  of  motion  of  the  atoms  for 
the  14th  normal  mode. 

(4)  The  geometry  of  the  system  can  be  changed  in  the  DRAW  program 
using  the  EDIT  command.  In  the  EDIT  mode  atoms  can  be  removed  or  changed  into 
atoms  of  other  elements.  Bond  lengths,  angles,  and  dihedrals  can  be  changed, 
and  atoms  can  be  added.  The  new  geometry  is  drawn  immediately  upon  QUITting 
the  EDIT  mode.  A  data  file  corresponding  to  the  new  geometry  is  generated  by 
issuing  the  OUTPUT  command. 
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We  have  now  seen  how  the  DRAW  program  can  be  used  to  help  interpret  and 
present  MOPAC  results.  In  the  next  section  we  will  work  through  some  examples 
of  MOPAC  calculations  that  will  help  to  improve  your  understanding  of  the 
capabilities  we  have  discussed  so  far.  We  will  also  introduce  some  new 
functions  and  approaches  that  should  help  you  begin  to  appreciate  the 
wide-ranging  powers  of  this  program. 

5.  Examples  of  applications. 

a.  Napthalene  and  the  keyword,  SYMMETRY. 

The  data  file  and  the  DRAW-generated  geometry  of  this  file  for  naphthalene 
are  shown  in  Fig.  3.  You  should  satisfy  yourself  that  the  geometry  shown  is 
consistent  with  the  geometry  definition  in  the  data  file.  The  keyword, 
SYMMETRY,  and  the  additional  set  of  numbers  below  the  geometry-ending  "blank 
line"  are  the  only  new  features  beyond  what  has  already  been  Introduced. 

Our  knowledge  of  molecular  structure  tells  us  that  many  of  the  geometric 
features  of  this  molecule  are  equivalent  due  to  its  symmetry.  For  example,  we 
know  that  the  bonds  3-2,  7-2,  5-1,  and  8-1  must  all  be  the  same  length  in  the 
equilibrium  geometry.  We  also  know  that  the  angles  3-2-1,  5-1-2,  7-2-1,  and 
8-1-2  must  be  equal.  For  reasons  which  we  will  discuss  later,  it  is 
advantageous  to  use  this  knowledge  by  optimizing  one  of  the  coordinates  in 
each  of  these  equivalent  sets  and  fixing  the  others  in  the  set  at  the 
optimized  value. 

In  this  example  the  bond  length  3-2  is  flagged  with  a  "1",  indicating  that 
this  coordinate  is  to  be  optimized.  The  distances  7-2,  5-1,  and  8-1  are 
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t  =3600  SYMMETRY 
naphthalene 


0.0  0 

1.4  1 

1.4  1 

1.4  1 

1.4  0 

1.4  0 
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0  0 
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.0  0 


Fig.  3.  Naphthalene  geometry  definition  with  symmetry. 
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flagged  with  "0"  to  indicate  that  they  are  not  optimized  by  the  normal 
procedure.  Instead,  we  will  include  symmetry  definitions  in  the  data  file  to 
assign  them  the  same  value  as  the  3-2  optimized  value.  How  is  this  done?  Note 
the  first  three  lines  of  the  symmetry  numbers  below  the  "blank  line."  The 
first,  "3  I  5",  tells  MOPAC  that  the  distance  coordinate  for  atom  5  (in  this 
example,  5-1,  as  seen  from  the  connectivity  list)  must  be  assigned  the  same 
value  as  the  optimized  value  for  the  distance  coordinate  for  atom  3  (in  this 
example,  3-2,  as  seen  from  the  connectivity  list).  The  number  "1"  in  between 
the  "3"  and  the  "5"  on  the  line  is  a  code  that  identifies  the  equivalent 
coordinates  as  a  distance.  The  next  two  lines,  "3  1  7"  and  "3  1  8",  say  the 
same  rule  applies  to  the  7-2  and  8-1  distances.  The  seventh  through  ninth 
lines  In  the  symmetry  numbers,  "3  2  5",  "3  2  7",  and  "3  2  8",  have  the  code 
number  "2"  for  "angle"  in  the  middle.  Thus  the  angle  3-2-1  will  be  optimized 
(as  Indicated  by  the  "1"  flag  in  the  geometry  definition)  and  the  angles 
5-1-2,  7-2-1,  and  8-1-2  (all  of  which  are  flagged  with  "0")  will  be  assigned 
the  optimized  3-2-1  value. 

To  summarize  what  we  did  in  this  specific  example  of  applying  the  SYMMETRY 
function:  the  middle  number  of  each  line  indicates  which  type  of  coordinate 
is  being  "symmetrized,"  the  first  number  tells  which  atom  will  have  that 
coordinate  optimized,  and  the  third  number  tells  which  atom  will  have  that 
coordinate  set  equal  to  the  equivalent  coordinate  of  the  optimized  atom.  In 
other  applications,  the  atom  corresponding  to  the  first  number  (called  the 
"reference  atom")  may  not  be  optimized,  and  the  middle  number  may  represent  a 
more  complex  relation  than  simple  equality.  A  complete  description  of  the 
ways  the  SYMMETRY  functions  can  be  used  is  given  in  Chapter  2  of  the  MOPAC 
Manual . 
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Note  that  the  numbering  of  the  atoms  in  this  example  is  not  the 
conventional  IUPAC  numbering  which  might  seem  more  desirable.  The  numbering 
scheme  was  chosen  to  allow  the  maximum  number  of  equivalent  sets  of 
coordinates  to  be  "symmetrized."  An  ability  to  find  the  best  numbering  scheme 
will  require  some  practice  with  the  SYMMETRY  function. 

But  why  use  symmetry?  First  it  can  save  significantly  on  the  time 
required  for  a  calculation,  a  factor  which  obviously  becomes  more  important 
with  larger  systems.  A  not  so  obvious  advantage  is  that  the  appearance  of 
exactly  equal  coordinates  in  the  tabulated  output  geometries  makes  the 
geometries  easier  to  read.  It  also  results  in  marginally  better  geometries 
and  heats  of  formation.  When  IR  spectra  are  calculated  from  "symmetrized" 
geometries,  degenerate  frequencies  come  out  to  be  almost  exactly  equal, 
whereas  differences  of  up  to  20  wavenumbers  may  result  when  "unsymmetrized" 
geometries  are  used.  More  advanced  uses  of  SYMMETRY  involve  finding 
transition  states  which  are  known  or  suspected  to  be  symmetrical  and  in 
defining  reaction  paths. 

b.  Vibrational  modes  and  IR  spectrum  of  A 1 2CI 7- - 

MOPAC  generates  many  results  which  can  be  compared  with  experimental 
data.  One  of  the  most  interesting  of  these  is  the  vibrational  frequencies  of 
a  system  calculated  by  using  the  keyword,  FORCE.  In  this  example  we  will 
demonstrate  how  this  function  Is  used  to  generate  the  IR  spectrum  of  the 
A 1 2C 1 7 ~  Ion. 
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We  begin  by  calculating  an  optimized  geometry  for  A^Cly  in  the  same  way 


we  did  the  CH3OH  molecule.  We  can  then  create  a  data  file  to  do  the  FORCE 
calculation  by  editing  the  resulting  .ARC  file.  We  simply  delete  all  of  the 
file  except  the  keywords  and  geometry  definition  which  appears  in  the  .ARC 
file  under  the  "FINAL  GEOMETRY."  We  add  the  keyword,  FORCE,  to  produce  the 
data  file  shown  below: 


FORCE  CHARSE— 1  T-1800 
CALC  OF  AL2CL7-  -  FORCE  CALC 


A1 

0.000000 

0 

0.000000 

0 

0.000000 

0 

0 

0 

0 

0.9403 

Cl 

2. 109916 

1 

0.000000 

0 

0.000000 

0 

1 

0 

0 

-0. 4071 

Cl 

2.118764 

1 

113.297154 

1 

0.000000 

0 

1 

2 

0 

-0.4220 

Cl 

2.118761 

1 

113.288224 

1 

129.847021 

1 

1 

2 

3 

-0.4218 

Cl 

2.299359 

1 

109.759524 

1 

-115.000056 

1 

1 

2 

3 

-0.3813 

A1 

2.297374 

1 

125.471304 

1 

-1.177332 

1 

5 

1 

2 

0.9419 

Cl 

2. 121224 

1 

101.504516 

1 

-179.671628 

1 

6 

5 

1 

-0.4249 

Cl 

2.112476 

1 

107.241245 

1 

-61.477832 

1 

6 

5 

1 

-0.4125 

Cl 

2.112748 

1 

107. 168291 

1 

62. 123575 

1 

6 

5 

1 

-0.4126 

0 

0.000000 

0 

0.000000 

« 

0 

0 

0 

0 

(NOTE:  The  charges  on  the  atoms  printed  in  the  extreme  right  column,  which 
come  with  the  .ARC  file,  do  not  need  to  be  deleted  for  the  data  file  to  work.) 


Performing  a  MOPAC  calculation  on  this  file  results  in  an  output  file 
which  initially  presents  the  same  kind  of  information  we  discussed  for  the 
CH30H  output  file.  After  the  INTERATOMIC  DISTANCES,  we  see  a  new  listing,  the 
INTERNAL  COORDINATE  DERIVATIVES,  which  shows  how  the  energy  varies  with 
changes  in  the  internal  coordinates  of  the  optimized  geometry.  These  should 
(and  almost  certainly  will)  be  less  than  0.5  kcal/mol/A.  Closely  following 
this  are  two  listings  of  the  "PRINCIPAL  MOMENTS  OF  INERTIA"  and  the 
"ORIENTATION  OF  MOLECULE  IN  FORCE  CALCULATION."  The  latter  gives  the 
cartesian  coordinates  of  the  atoms  with  the  system  center  of  gravity  at  the 
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zero  of  the  coordinate  system.  Next  comes  a  list  of  the  steps  used  to  set  up 
the  Hessian  which  Is  diagonalized  to  get  the  force  constants.  A  condensed 
form  of  the  Hessian  follows  In  the  table  designated  "FORCE  MATRIX  IN 
MILLIDYNES/ANGSTROMS . "  The  numbers  In  this  table  Indicate  the  strength  or 
"stiffness"  of  the  bonds  connecting  the  atoms  in  the  molecule.  Next  is  a 
listing  of  the  results  of  the  "NORMAL  COORDINATE  ANALYSIS."  This  is  of 
limited  value  for  most  users. 


Finally  we  get  to  the  key  result,  the  "DESCRIPTION  OF  VIBRATIONS."  The 
first  few  vibrations  may  not  be  listed  if  their  frequencies  are  less  than  50 
cm-1  since  these  tend  to  have  large  uncertainties.  Those  modes  with  higher 
frequencies  are  listed  in  order  of  increasing  frequency.  We  will  discuss  only 
the  main  features  of  this  listing;  others  are  discussed  in  Chapter  5  of  the 
MOPAC  Manual . 


VIBRATION 

14 

ATOM  PAIR 

« NEWS'/  CONTRIBUTION 

RADIAL 

FREQ. 

313.  80 

Cl  3  —  A1  6 

20. IX  (  99.6X) 

64.2X 

T-DIPOLE 

1 . 0483 

A1  1  —  Cl  5 

19.57. 

61.97. 

TRAVEL 

0.0800 

A1  1  —  Cl  2 

11.07. 

33.  9X 

RED.  MASS 

16. 7786 

A1  6  —  Cl  8 

10.  7X 

28. 6X 

A1  6  —  Cl  9 

10.77. 

28.  3X 

VIBRATION 

IS 

ATOM  PAIR 

ENERSv  CONTRIBUTION 

RADIAL 

FREQ. 

348. 16 

Al  1  —  Cl  3 

33. 9X  (  78.97.) 

86.67. 

T-DIPOLE 

4.9923 

Cl  3  —  Al  6 

35.  77. 

86.  3X 

TRAVEL 

0.  0793 

Al  1  --  Cl  2 

4.9X 

17.  0X 

RED.  MASS 

17.0960 

Al  1  —  Cl  4 

4.BX 

20.87. 

The  "FREQ"  In  cm-1  and  "T-DIPOLE"  (transition  dipole  in  Debyes)  are  given 
in  the  left  column  as  seen  for  the  examples  shown  above  for  the  vibrations,  14 
and  15,  for  AI2CI7-.  The  next  column  lists  all  pairs  of  atoms  that  contribute 
significantly  to  the  energy  of  the  mode.  In  the  next  column  is  an  indication 
of  how  much  each  pair  contributes  to  the  energy  of  the  mode.  A  discussion  of 
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how  these  numbers  are  calculated  and  what  they  strictly  mean  Is  given  In  the 
Manual.  For  the  beginning  user,  the  following  may  suffice.  The  energy  of  a 
photon  stimulating  a  given  mode  Is  distributed  throughout  the  molecule.  Pairs 
of  atoms  with  the  highest  "ENERGY  CONTRIBUTION"  are  seen  as  absorbing  most  of 
that  photon's  energy.  The  number  In  parentheses  following  the  top  pair  Is  the 
energy  of  that  pair  divided  by  the  total  energy  of  the  mode.  If  the  78. 9X 
given  for  the  A11-C15  pair  in  vibration  15  seems  high  (especially  since  the 
second  pair,  C15-A16,  has  about  the  same  value  for  the  ENERGY  DISTRIBUTION), 
it  need  only  be  recognized  that  both  pairs  share  a  common  atom,  C15,  so  that 
in  effect  there  is  double  counting  in  this  method  for  describing  the  way  the 
energy  Is  distributed.  (It  is  probably  small  consolation  that  other  methods 
of  description  are  even  less  satisfying!) 

The  last  column  gives  a  measure  of  how  much  the  motion  of  the  atoms  In  the 
pair  Is  radial  (as  occurs  in  pure  stretching)  for  the  mode.  The  smaller  this 
number  is,  the  more  the  motion  between  the  atoms  can  be  characterized  as 
bending  or  wagging.  This  information  is  useful  in  describing  the  mode  in 
traditional  terms. 

Vibrations  14  and  15  can  be  reasonably  assumed  to  be  the  symmetric  and 
asymmetric  stretches  in  the  A1-C1-A1  bridge.  Diagrams  of  these  modes  can  be 
obtained  by  using  DRAW  on  the  output  (not  .ARC)  file.  The  command,  DISPLAY 
VIBRATION  14  QUIT,  in  DRAW  gives  the  picture  of  the  mode  shown  in  Fig.  4. a. 
Vibration  15  is  shown  in  Fig.  4.b.  Fig.  5  shows  the  calculated  frequencies 
and  their  intensities  (obtained  as  T-DIPOLE  squared  multiplied  by  the 
frequency)  for  AI2CI7"  compared  to  the  IR  spectrum  obtained  experimentally. 
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Fig.  4.  A 1 2CI 7“  vibrations 
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iinteni 


*  Hvistendahl ,  J.;  Klaeboe,  P.;  Rytter,  E.;  0ye,  H.  A.  Inorg.  Chem.  1984 ,  23,  p.  706. 

Fig.  5.  A 1 2C 1 7 -  calculated  vs  experimental  IR  spectra. 
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c.  Reaction  path  calculation. 

In  this  example  we  consider  the  problem  of  finding  an  energy  of  activation 
for  the  reaction, 

Br-  +  CH3CI  — >  CH3Br  +  Cl"  . 

The  key  to  determining  an  activation  energy  is  identifying  and  calculating 
the  AHf  for  the  transition  state  for  the  reaction.  Usually  in  MOPAC 
calculations,  this  involves  a  three  step  process.  First,  we  get  an 
approximate  geometry  for  the  transition  state  by  finding  the  maximum  in  a 
reaction  profile  (i.e.,  a  plot  of  AHf  for  the  reacting  system  vs.  the  selected 
reaction  coordinate).  Next  we  refine  that  geometry  by  locating  a  stationary 
point  at  that  maximum  (I.e.,  a  point  at  which  the  first  derivative  of  the 
energy  with  respect  to  the  coordinates  approaches  zero).  Finally  we 
characterize  that  stationary  point  to  determine  if  it  is  indeed  a  transition 
state.  We  will  now  briefly  describe  how  each  of  these  step'  esn  be  carried 
out  for  our  example. 

(1)  Finding  maximum  in  reaction  profile. 

We  begin  by  creating  a  data  file  for  the  initial  geometry  of  the  reacting 
system  as  shown  below. 

T-72B0  CHARGE— 1 

REAX  OF  BR-  WITH  CH3CL 

C  0  0  0  0  0  0  0  0  0 

cl  1 .  e  1  0  0  00  100 

H  1.2  1  109.0  1  0  0  1  2  0 

H  1.2  1  109.0  1  120.0  1  123 

H  1.2  1  109.0  1  -120.0  1  123 

BR  20.0  -1  180.0  1  0.0  1  123 

10.0  3.0  4.0  3.3  3.0  2.9  2.8  2.7  2.6  2.5  2.4  2.3  2.2  2.1  2.0  1.9  1.8 
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The  keyword,  CHARGE=-1,  gives  the  total  charge  of  the  system,  which  initially 
(with  the  Br~  well  removed  from  the  CH3CD  resides  entirely  on  the  Br~.  All 
the  coordinates  in  the  geometry  description  are  flagged  with  "1",  indicating 
they  are  to  be  optimized,  except  for  one.  This  one  coordinate,  the  Br~-C 
distance,  has  been  chosen  to  be  the  "reaction  coordinate."  It  is  flagged  with 
a  "-1"  to  indicate  that  it  will  be  initially  fixed  at  the  20.0  A  specified  in 
the  geometry  definition.  AHf  for  this  initial  configuration  will  be 
calculated  for  the  geometry  in  which  all  the  other  coordinates  have  been 
optimized.  Then  the  Br~  will  be  fixed  at  shorter  distances  from  the  C  and  the 
process  will  be  repeated.  The  values  at  which  the  Br~-C  distance  will  be 
fixed  are  specified  in  the  line  immediately  below  the  "blank  line"  at  the  end 
of  the  initial  geometry  definition.  The  Br--C  distances  should  be  decreased 
to  beyond  the  value  it  would  be  expected  to  have  when  the  reaction  has 
occurred.  The  intent  here  is  to  generate  a  reaction  profile  (that  is,  a  plot 
of  the  system  heat  of  formation  as  a  function  of  the  reaction  coordinate) 
which  passes  through  a  maximum  (see  Fig.  6. a.).  The  transition  state  is  then 
assumed  to  be  near  that  maximum  and  further  efforts  can  be  made  to  locate  it 
and  characterize  it  as  a  transition  state. 

The  output  from  this  calculation  is  essentially  a  sequence  of  sets  of 
results  similar  to  the  output  from  the  calculation  we  did  previously  on 
CH3OH.  At  the  very  beginning  of  the  output,  the  "POINTS  ON  REACTION 
COORDINATE"  are  listed.  These  are  the  values  we  selected  to  have  the  Br_-C 
distance  fixed  at.  There  will  be  a  set  of  results  for  ea^h  fixed  Br~-C 
distance,  and  these  can  be  interpreted  in  the  same  way  as  previously  discussed 
for  the  CH30H  output. 


If  the  AHf  values  are  plotted  as  a  function  of  these  fixed  distances  along 
the  "reaction  coordinate,"  we  obtain  the  curve  shown  in  Fig.  6(a).  Note  that 
there  is  a  minimum  in  the  curve  before  the  maximum  is  reached.  This 
corresponds  to  the  Van  der  Waals  complex  formed  by  these  two  species.  The 
three  darkened  points  in  the  vicinity  of  the  maximum  indicate  where  we  must 
look  to  find  the  transition  state.  Thus  we  next  perform  another  reaction  path 
calculation  (with  smaller  steps)  starting  with  the  first  of  these  points,  at 
2.4  A,  and  ending  with  the  last  point  at  2.2  A.  We  can  create  the  data  file 
by  editing  the  .ARC  file  from  the  first  reaction  path  calculation.  We  simply 
delete  from  this  file  everything  except  the  keywords  and  geometry  definition 
of  the  "FINAL  GEOMETRY"  obtained  for  the  2.4  A  calculation.  Next  we  add, 
after  the  "blank  line,"  the  points  on  the  reaction  path  going  through  the 
maximum  as  shown  in  the  data  file  below: 

T- 7200  CHARGE— 1 
REAX  OF  BR-  WITH  CH3CL 

c  a.  mama  a  a.aaaaaa  a  a.  mama  a  a  a  a  0.3632 

Cl  1 . 972184  1  0.000000  0  0.000000  0100  -0.5739 

H  1.094351  1  97. 129B43  1  0.000000  0120  0.0233 

H  1.094488  1  97.227117  1  120.004835  1123  0.0234 

H  1.094510  t  97.213600  1  -120.033132  1123  0.0234 

Br  2.400000  -1  179.917373  1  -5.61490B  1123  -0.B393 

2.38  2.36  2.34  2.32  2.30  2.28  2.26  2.24  2.22  2.20 

The  results  of  this  calculation  are  plotted  in  Fig.  6(b).  We  are  now  In  a 
position  to  make  an  attempt  to  confirm  what  appears  to  be  a  valid  transition 
state . 

(2)  Locating  a  stationary  point  at  the  maximum. 

We  choose  the  point  just  before  the  maximum  in  the  curve,  at  2.26  A,  as 
the  geometry  which  we  will  now  allow  to  fully  optimize.  This  means  we  will 
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Fig.  6.  Br-  +  CH3CI  reaction  profiles. 
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not  fix  the  8r~-C  distance.  However  we  must  also  tell  MOPAC  that  we  do  not 


want  to  find  the  minimum  in  &Hf  from  this  starting  geometry,  which  would  take 
us  back  down  the  reaction  profile.  What  we  want  is  the  maximum  on  the 
reaction  profile,  i.e.,  the  transition  state.  To  do  this,  we  must  use  a 
keyword  <we  will  use  NLLSQ)  which  tells  MOPAC  to  minimize  the  GNORM,  the 
derivative  of  AHf  with  respect  to  the  coordinates.  (Other  keywords  which  can 
be  used  to  locate  the  transition  state  are  SIGMA  and  SADDLE,  and  these  are 
discussed  in  the  MOPAC  Manual.)  Again  we  edit  the  .ARC  file  from  this  last 
calculation  to  create  a  data  file  which  has  the  keywords  and  geometry 
definition  of  interest;  in  this  case,  that  for  the  2.26  A  point.  We  change 
the  "-1"  flag  in  the  .ARC  file  to  a  "1"  for  the  Br--C  distance  to  indicate 
that  it  will  be  optimized  along  with  the  other  coordinates.  Then  we  add  the 
keyword,  NLLSQ.  If  the  calculation  is  successfully  concluded,  the  resulting 
GNORM-optimized  geometry  and  AHf  correspond  almost  certainly  to  the  transition 
state. 


(3)  Characterization  of  the  stationary  point. 

We  now  test  whether  the  stationary  point  obtained  in  Step  (2)  is  indeed  a 
transition  state.  The  final  geometry  is  extracted  out  of  the  .ARC  file  and 
another  data  file  created.  The  NLLSQ  keyword  should  be  removed  and  the 
keyword,  FORCE,  used.  This  gives  the  frequencies  of  the  vibrations  for  the 
stationary  state.  If  the  transition  state  is  valid,  there  will  be  one,  and 
only  one,  vibrational  mode  with  a  negative  frequency.  In  our  example, 
vibration  1  is  seen  to  be  the  only  mode  with  a  negative  frequency.  Fig.  7 
shows  a  DRAW-generated  picture  of  this  mode,  which  effectively  describes  the 
actual  reaction  coordinate  as  well  as  the  geometry  of  the  transition  state. 
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Br~  +  CH3CI  transition  state. 
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Note  that  the  vibration  corresponds  to  the  dissociation  of  the  transition 
state  into  Ci^Br  and  Cl". 

Once  the  transition  state  has  been  fully  confirmed,  the  energy  of 
activation  can  be  obtained  as  the  difference  in  AHf  between  the  fully 
separated  species  and  the  transition  state.  It  is  left  as  a  useful  exercise 
to  find  the  transition  state  for  the  reverse  reaction, 

Cl"  +  CH3Br  —  >  CH3C1  +  Br"  . 

6.  Other  considerations. 

We  now  discuss  some  other  aspects  of  the  MOPAC  program  which  are  not 
essential  to  operation  at  the  Introductory  level  but  are  very  useful  even  to 
the  beginning  user. 

a.  Restart  files. 

If  a  job  uses  more  than  one  hour  of  CPU  time,  a  restart  file,  called 
FILENAME. RES,  Is  created  which  saves  the  results  which  have  been  calculated  in 
that  hour.  A  file  containing  the  density  matrix,  FILENAME. DEN,  is  also 
created.  After  each  additional  hour,  a  new  set  of  these  files  is  created  to 
replace  the  previous  one.  If  for  any  reason  the  calculation  is  terminated 
(e.g.,  a  computer  shutdown,  failure  to  achieve  a  self-consistent  field),  the 
calculation  can  be  picked  up  from  the  point  at  which  the  last  restart  file  was 
created.  This  can  be  accomplished  by  adding  the  keyword,  RESTART,  to  the  data 
file  for  the  calculation  and  resubmitting  it.  The  geometry  definition  must 
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not  be  changed  In  this  resubmitted  file.  If  the  calculation  was  stopped 
because  of  some  problem  with  the  calculation  Itself,  the  problem  would 
normally  be  described  at  the  end  of  the  output  file.  In  this  case  the 
keywords,  RESTART  and  1SCF  (which  tells  MOPAC  to  perform  only  one  SCF  and  then 
stop),  should  be  added  to  the  data  file.  The  .ARC  file  output  from  this 
calculation  gives  the  geometry  before  the  problem  with  the  calculation  forced 
it  to  end.  This  file  probably  can  be  edited  to  create  a  new  data  file  which 
avoids  the  calculational  problem. 

b .  Dummy  atoms . 

A  dummy  atom  is  a  pure  mathematical  point  which  can  be  entered  into  a 
geometry  definition  in  the  same  way  that  other  atoms  are  entered.  Coordinates 
for  other  atoms  can  then  be  specified  by  reference  to  the  dummy  atom.  The 
symbol  for  a  dummy  atom  is  XX.  Dummy  atoms  are  used  when  their  presence  as  a 
point  of  reference  simplifies  a  geometry  definition.  This  may  be  especially 
useful  when  there  are  more  than  one  species  in  the  system  being  calculated. 

It  can  be  almost  indispenslble  in  the  calculation  of  specific  reaction  paths. 
For  example,  suppose  it  is  desired  to  follow  the  reaction  of  Cl 2  and  CH2=CH2 
with  the  C=C  and  Cl 2  bonds  approaching  each  other  while  remaining  parallel  to 
each  other.  Dummy  atoms  centered  in  the  C=C  and  Cl 2  bonds  greatly  simplify 
this  reaction  path  calculation. 

c.  Sparkles. 

Sparkles  are  computational  constructs  made  to  simulate  unpolarizable  ions 
of  diameter  1.4  A.  The  charge  of  the  sparkle  is  set  by  the  symbol  used  to 
enter  it  into  the  geometry  definition:  +,  ++,  -,  or  — .  They  have  no  heats 
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of  atomization,  no  ionization  potentials,  and  no  orbitals.  They  are 
computationally  designed  so  that  other  atoms  or  sparkles  can  approach  them  to 
within  about  2  A  quite  easily  but  can  come  closer  than  1.4  A  only  with  great 
difficulty.  The  most  obvious  use  of  sparkles  is  to  simulate  counterions. 

Other  uses  are  discussed  In  Chapter  6  of  the  Manual. 

d.  Condensed  phase  systems. 

The  species  in  a  MOPAC  model  are  surrounded  by  vacuum  so  that  calculated 
results  most  closely  approximate  the  behavior  of  the  "real"  species  in  the  gas 
phase.  MOPAC  results  can  be  applied  to  liquid  phase  system,  but  only  with 
great  caution.  In  aqueous  systems,  solvated  species  will  in  most  cases  have 
significantly  different  properties  than  what  they  have  in  the  gas  phase.  This 
can  be  attributed  to  the  high  dielectric  constant  of  water  which  significantly 
modifies  electrostatic  interactions.  Furthermore,  solute-solvent  interactions 
change  certain  properties  (or  even  the  chemical  nature)  of  both  solute  and 
solvent  species.  As  a  fairly  obvious  example,  MOPAC  predicts  the  system,  NH3 
+  HCOOH ,  to  be  stable  In  the  gas  phase.  However,  in  aqueous  solution,  we  know 
that  solvated  NH4+  and  HCOO-  are  the  stable  species.  In  some  cases,  addition 
of  solvent  molecules  to  the  model  can  be  tried  in  order  to  predict  behavior  in 
solution,  but  this  must  be  done  very  carefully.  For  species  in  low  dielectric 
constant  solvents,  the  effect  on  accuracy  of  not  including  solvent  species  In 
the  model  is  less.  Some  properties,  such  as  certain  vibrational  frequencies, 
can  be  fairly  insensitive  to  phase,  as  experimental  data  have  shown. 
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e.  A  final  word. 

The  best  way  to  develop  your  ability  to  use  MOPAC  effectively  is  to  use  it 
on  systems  with  which  you  are  familiar.  When  problems  with  interpretation  or 
apparent  bugs  in  the  program  arise  which  you  cannot  solve  using  the  MOPAC 
Manual,  you  are  encouraged  to  contact  Dr.  Jimmy  Stewart  at  FJSRL. 

Write  to:  Dr.  J.  J.  P.  Stewart 

FJSRL/NC 

US  Air  Force  Academy 

Colorado  Springs  CO  80840-6528 

USA 

or  call  at:  (303)  472-2655  or  AV  259-2655. 

or:  stewart@usafa.arpa  on  the  arpa  network. 

MOPAC  (the  Manual  and  program)  Is  revised  on  a  nearly  annual  basis  and  the 
current  version  is  always  available  from  the  Quantum  Chemistry  Program 
Exchange  (QCPE);  Department  of  Chemistry;  Indiana  University;  Bloomington, 
Indiana  47405.  Contact  the  Editor,  Richard  Counts  at  (812)  335-4784  for 
details. 
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APPENDIX  A 


More  on  Geometry  Definition 

The  day  will  come  when  a  chemist  can  construct  a  ball  and  stick  model  for 
the  molecule  he  wants  to  perform  calculations  on,  put  the  model  into  an  image 
processor,  and  get  out  a  geometry  definition  which  his  software  can  work  with 
as  a  starting  geometry.  Actually  certain  stereoscopic  graphic  systems  already 
allow  you  to  draw  3-D  pictures  which  can  be  converted  to  data-files,  but  these 
are  not  yet  completely  routine  lab  tools.  Until  then,  it  is  necessary  to  know 
how  to  construct  a  geometry  definition  in  internal  coordinates  from  the  ball 
ind  stick  model.  Here  we  amplify  somewhat  the  description  of  how  to  construct 
a  geometry  definition  which  is  given  in  Section  3. a. 

Fig.  A-l  is  a  graphical  representation  of  the  process  by  which  the  three 
internal  coordinates  needed  to  specify  an  atom,  I,  are  determined.  "Distance" 
and  "angle"  are  readily  understood  once  it  is  clear  that  I  represents  the  atom 
whose  coordinates  are  being  specified  and  that  A,  B,  and  C  refer  to  previously 
defined  atoms  in  the  molecule.  When  taking  guesses  at  initial  distances  where 
I  and  A  are  actually  bonded  together,  it  is  generally  better  to  slightly 
overestimate  the  bond  length. 

The  barrel-and-sight  method  is  recommended  as  an  easily  learned  way  to 
determining  what  the  value  and  sign  of  the  "dihedral"  angle  should  be.  In 
Fig.  A-l,  the  dihedral  is  shown  as  the  angle  between  the  I-A-B  and  A-B-C 
planes.  The  intersection  of  the  two  planes,  i.e.  B-A,  can  be  considered  a 
barrel  down  which  the  viewer  sights  in  order  to  determine  the  dihedral.  The 
sight  is  the  line  which  joins  atom  C  to  the  B-A  "barrel.”  With  the  eye 
looking  down  the  barrel,  the  position  of  I  relative  to  the  sight  can  be 
evaluated  in  terms  of  how  much  of  an  angle  it  is  rotated  away  from  the  sight 
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and  whether  the  rotation  is  clockwise  (CW)  or  counterclockwise  (CCW).  If  it's 
CW,  the  dihedral  is  positive;  if  it's  CCW,  the  dihedral  is  negative.  Thus 
instead  of  entering  a  dihedral  of  270  degrees,  one  would  enter  it  as  -90 
degrees.  As  shown  in  Fig.  A-l,  the  eye  msut  be  at  the  B  end  of  the  barrel. 

It  must  also  be  kept  in  mind  that  the  atoms,  I,  A,  B,  and  C,  in  a  dihedral 

definition  do  not  have  to  be  bonded  to  each  other.  Finally  it  must  be 
remembered  that  the  dihedral  I-A-B-C  cannot  be  defined  if  atoms  A,  B,  and  C 
are  in  a  straight  line.  Occasionally,  this  condition  occurs  during  a  geometry 
optimisation  and  some  change  in  geometry  definition  must  be  made  to  avoid  this 
"job-stopper. " 

In  deciding  on  the  order  in  which  the  atoms  in  the  geometry  are  defined, 
it  should  be  kept  in  mind  that  atom  1  will  be  the  origin  of  the  system  in 

cartesian  coordinates  and  atom  2  will  be  along  the  x-axis,  with  its  y  and  z 

coordinates  at  zero.  Thus  the  I-A  =  1-2  distance  will  be  a  horizontal  in  the 
plane  of  the  screen  when  displayed  through  the  DRAW  program.  Although  you  can 
use  DRAW  commands  to  orientate  the  molecule  any  way  you  want  from  this 
starting  position,  it  is  convenient  to  have  a  desired  orientation  appear 
immediately  by  using  an  appropriate  ordering  of  atoms  in  the  geometry 


definition. 


SPECIFY 


BY  REFERENCE  TO 


Distance  : 
Angle  : 

Dihedral  : 


I  -  A  - 


B 

B  -  C 


IN  CONNECTIVITY  LIST 


BARREL 

Figure  A-1 
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APPENDIX  B 


MOPAC  Calculations  on  Polymers 


While  useful  results  can  be  obtained  for  polymers,  it  must  be  remembered 
that  in  MOPAC,  polymers  are  treated  as  single,  untangled,  non-interacting 
molecules  of  infinite  length.  If  you  know  how  to  create  a  data  file  for  a 
molecule,  you  can  create  a  data  file  for  a  polymer  by  using  just  two 
additional  concepts:  the  translation  vector  (TV)  and  the  polymer  unit  cell 


(PUC). 

1.  The  PUC  is  a  segment  of  the  polymer  which  will  be  included  in  the 
geometry  definition  of  the  data  file.  It  must  contain  at  least  the  complete 
repeating  unit  (let's  call  this  the  monomer).  However,  it  may  contain  more 
than  one  monomer  because  it  must  be  at  least  8  A  in  length.  A  valid  PUC  for 
polyethylene  is  shown  below: 


TRANSLATION  VECTOR 


TV(26) 

i 

POSITION 
OF  XX(9) 


Note  that  the  PUC  contains  4  monomers.  The  geometry  definition  will  also  have 
a  dummy  atom,  XX,  in  the  position  that  the  carbon,  C(l),  would  occupy  if  the 
PUC  were  to  be  repeated.  Thus  XX  is  at  the  equivalent  point  of  CO)  in  the 
adjacent  PUC. 

2.  The  TV  is  a  distance  vector  which  "translates"  a  point  in  the  PUC  to 
the  equivalent  point  in  the  adjacent  PUC.  Thus  the  drawing  above  shows  the  TV 
which  takes  CO)  to  its  equivalent  point  (at  which  XX  is  located)  in  the 
adjacent  PUC.  The  TV  must  be  greater  than  8  A  to  insure  that  atoms  in  the 
middle  of  the  PUC  do  not  Interact  simultaneously  with  atoms  at  both  ends 
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of  the  PUC.  This  Insures  that  all  atoms  In  the  PUC  behave  as  If  they  were 
well  removed  from  the  ends  of  the  polymer.  The  TV  Is  entered  as  the  distance 
coordinate  on  the  last  line  (before  the  blank  one)  of  the  geometry  definition. 

POLYETHYLENE  DATA  FILE.  The  data  file  for  the  polyethylene  model 
discussed  above  is  shown  on  the  attached  page.  Note  that  it  was  convenient  to 
enter  the  XX  (9)  equivalent  point  after  the  carbons  had  been  entered.  The 
geometry  definition  for  the  H's  is  entered  next.  The  value  for  the  TV  is 
obtained  by  taking  the  data  file  (without  the  TV  entered)  into  the  DRAW 
program.  In  DRAW,  the  INFORMATION  function  is  used  to  obtain  the  value  of  the 
C(1)  -  XX(9)  distance.  This  distance  is  entered  on  the  last  line  of  the 
geometry  definition  which  starts  with  "TV"  as  the  atom  symbol.  It  Is  flagged 
with  a  "1"  to  indicate  it  is  to  be  optimized.  0's  are  entered  for  the  angle 
and  dihedral  values  and  their  flags.  The  connectivity  values  for  TV  are 
assigned  as  fol lows : 

A  -  The  number  of  the  atom  which  is  being  "translated"  into  the 
adjacent  PUC. 

B  -  The  number  of  the  XX  atom  which  is  positioned  at  the  eqivalent 
point  for  the  atom  being  tanslated. 

C  -  Any  atom  previously  defined. 

This  data  file  is  entered  with  no  other  special  keywords  (the  keyword, 

XYZ ,  must  not  be  used).  The  attached  output  file  contains  some  explanatory 
remarks  on  certain  features  that  differ  from  other  output  files.  However,  as 
can  be  seen,  the  contents  do  not  differ  significantly  from  the  output  file  for 


a  normal  molecule. 


The  ways  In  which  MOPAC  calculations  on  polymers  can  be  used  are  still  being 
discovered.  Some  uses  currently  being  made  of  these  results  Include 
calculating  heats  of  polymerization, *  IR  Spectra,  and  elastic  moduli.** 

*  Stewart,  J.J.P.,  New  Polymeric  Mat.  1  (1)  53-61  (1987). 

**  Kiel,  H.E.,  and  Stewart  J.J.P.,  Int.  3.  Quantum  Chem;  Quantum  Chemistry 
Symposium  20,  529-540  (1986);  John  Wiley' and  Sons,  Inc. 
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DATA  FILE  FOR  POLYETHYLENE 


T-60000 

CALCULATION  OF  POLYETHYLENE 
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OPTIMIZED  GEOMETRY  OF  POLYETHYLENE  PUC 
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OUTPUT  FILE  FOR  POLYETHYLENE 

**********##***#******##*****#************#*****************************•******'* 
**  FRANK  J.  SEILER  RES.  LAB. t  U. S.  AIR  FORCE  ACADEMY,  COLO.  SP6S.  ,  CO.  80840  ** 
*#*********#**###***********##**********#***»********************************** 

MNDO  CALCULATION  RESULTS 
VERSION  4.0O 
Tue  Mar*  8  12:00:48  1988 

******************************************************************************* 

*  T*  -  A  TIME  OF  60000.0  SECONDS  REQUESTED 

*  DUMP*N  -  RESTART  FILE  WRITTEN  EVERY  3600.0  SECONDS 
************************************************************************************* 
T  *6000® 

CALCULATION  OF  POLYETHYLENE 


ATOM 
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BOND  ANGLE 

TWIST  ANGLE 
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£ 
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5 
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6 
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£4 
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* 
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8 

7 
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£5 

H 
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* 
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* 
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* 
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7 

6 

£6 

Tv 
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.  00000 
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9 

8 

THE  SYSTEM  IS  A  POLYMER 

Note  indication 

unit  cell  translation  vectors  of  polymer  calculation 

and  translation  vector. 

X  Y  Z  - 
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CARTESIAN  COORDINATES 
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.  0000 

.  0000 

.  0000 

£ 

C 

1. 5500 

.  0000 

.  0000 

3 

C 

2. 1804 

1. 4160 

.  0000 

A 

c 

3. 7304 

1. 4160 

.  0000 

c 

p 

4.  3609 

2. 8320 

|> 

u 

.  ic'Wvlci 

e 

7 

0 

c 

c 

c 

5.  9109 

6.  5413 

8. 0913 

2. 8320 

4.  £480 

4. 2480 

.  0000 
.  0000 
.  0000 

Note  that  XX  is  not 
included  in  calculation. 

9 

10 

1 1 

12 

H 

H 

H 

H  * 

-.3614 
-. 3614 

1. 8824 

1 . 8824 

.  5248 
.  5248 
5436 
-. 5436 

-.  9089 
.  9089 
.9089 
.  9089 

H  atoms 
with  first 
H(9). 

are  renumbered 
H  now  labelled 

13 

H 

1. 8480 

1. 9596 

.  9089 

14 

H 

1. 8480 

1. 9596 

-. 9089 

1 5 

H 

4. 0918 

.  8912 

9089 

16 

H 

4. 0918 

.  8912 

.  9089 

17 

H 

4. 0265 

3. 3756 

.  9089 

18 

H 

4. 0285 

3. 3756 

-. 9089 

19 

H 

6. 2723 

2. 3072 

9089 

£0 

H 

6.2723 

2. 3072 

.  9089 

£1 

H 

6. 2089 

4. 7916 

.  9089 

£2 

H 

6. 2089 

4. 7916 

9089 

£3 

H 

8. 4527 

3. 7232 

9089 

£4 

H 

8. 4527 

3. 7232 

.  9089 

H:  (MNDO)  :  M.  J. 

S.  DEWAR,  W.  THIEL,  J.  AM. 

CHEM. 

SOC 

.  ,  99, 

4899, 

(1977) 

C:  (MNDO):  M.  J. 

S.  DEWAR,  W.  THIEL,  J.  AM. 

CHEM.  SOC 

.  ,  99, 

4899, 

(1977) 

RHF  CALCULATION,  NO.  OF 

DOUBLY  OCCUPIED  LEVELS  =  £4 

INTERATOMIC 

DISTANCES 

NOT  SHOWN 

CYCLE: 

1 

TIME:  1 

34.  75 

TIME  LEFT: 

59815.  4 

GRAD.  : 

59. 898 

HEAT: 

-36. 15304 

CYCLE: 

£ 

TIME: 

94.  £5 

TIME  LEFT: 

59721 . 1 

GRAD. : 

16.  547 

HEAT  : 

-36. 86357 

CYCLE: 

3 

TIME: 

54.  45 

TIME  LEFT: 

59666. 6 

GRAD.  : 

13. 389 

HEAT: 

-37. 05841 

CYCLE : 

4 

TIME: 

87.  80 

TIME  LEFT: 

59578. 9 

GRAD.  : 

1 1 . 855 

HEAT: 

-37. 15977 

CYCLE : 

5 

TIME: 

54.  50 

TIME  LEFT: 

59524. 4 

GRAD.  : 

10.  413 

HEAT: 

-37. 28218 

CYCLE: 

6 

TIME: 

54.  45 

TIME  LEFT: 

59469. 9 

GRAD. : 

10.  176 

HEAT  : 

-37. 33411 

C'CLE : 

7 

TIME: 

87.  95 

TIME  LEFT: 

59382.  0 

GRAD.  : 

9.  275 

HEAT: 

-37. 39660 

CYCLE: 

0 

TIME: 

52.  95 

TIME  LEFT: 

59329.  0 

GRAD.  : 

8.  302 

HEAT  : 

-37. 40682 

CYCLE: 

9 

TIME:  102.90 

TIME  LEFT: 

59226.  1 

GRAD.  : 

6.  1  17 

HEAT: 

-37. 42818 

CYCLE: 

10 

TIME: 

52.  85 

TIME  LEFT: 

59173.  3 

GRAD.  : 

5.  044 

HEAT: 

-37. 43771 

CYCLE: 

1 1 

TIME:  1 

23.  35 

TIME  LEFT: 

59049.  9 

GRAD.  : 

3.  647 

HEAT : 

-37. 45544 

CYCLE : 

12 

TIME: 

54.  95 

TIME  LEFT: 

58995. 0 

GRAD.  : 

2.  70S 

HEAT: 

-37. 46193 

CYCLE : 

13 

TIME: 

84.  70 

TIME  LEFT: 

58910. 3 

GRAD.  : 

1 . 733 

HEAT: 

-37. 46963 

0TEST  ON 

X 

SATISFIED 

HOWEVER,  A  COMPONENT  OF 

GRADIENT  IS  LARGER  THAN 

.  12 

CYCLE: 

14 

TIME: 

49.  80 

TIME  LEFT: 

58860.  4 

GRAD.  : 

1 . 661 

HEAT  : 

-37. 47034 

HERBERTS  TEST  SATISFIED  -  GEOMETRY  OPTIMIZED 
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T -60000 

CftLCUL.PT  I  ON  OF  POLYETHYLENE 


HERBERTS  TEST  WAS  SATISFIED  IN  FLETCHER-POWELL 
SCF  FIELD  WAS  ACHIEVED 


MNDO 


CALCULATION 


VERSION 


This 

is  heat  of  formation 
of  the  PUC. 

T ue  Mar 

FINAL 

HEAT 

OF  FORMATION 

= 

-37. 47034 

kcal 

ELECTRONIC 

ENERGY 

-9996. 15906 

EV 

CORE-CORE  REPULSION 

= 

8744. 597£0 

EV 

IONIZATION 

POTENTIAL 

rr 

1  1.  55576 

4.  00 

a  12:  22-:  1 0  1988 


NO.  OF  FILLED  LEVELS 
MOLECULAR  WEIGHT 


SCF  CALCULATIONS 
COMPUTATION  TIM* 


£4 

1 1 £. £14 


71 


1 15£. 90  SECONDS 


Optimized  geometry. 
Shows  optimized 
positions  of 

XX  and  TV. 


ATOM 

CHEMICAL 

BOND  LENGTH 

BOND  ANGLE 

TWIST  ANGLE 

V 

lUMBER 

SYMBOL 

(ANGSTROMS) 

(DEGREES) 

(DEGREES) 

<I> 

NA:  I 

NB:NA: I 

NC:NB:NA: 

I 

NA 

NB 

1 

£ 

C 

C 

1.54345 

# 

1 

3 

C 

1.  54350 

* 

114. £0942 

■* 

2 

1 

4 

C 

1 . 54329 

* 

114. £0773 

* 

180.  00952 

* 

U 

2 

5 

C 

1. 543£9 

* 

114. £1747 

* 

179. 6481 1 

* 

4 

u 

6 

C 

1 . 54388 

* 

114.  £1079 

* 

179.  93238 

* 

5 

4 

7 

C 

1. 54334 

* 

114.  £181  1 

* 

179. 68444 

* 

6 

5 

8 

C 

1.54351 

* 

1 14.  19561 

* 

-179.  91091 

♦ 

7 

£ 

9 

XX 

1. 54530 

* 

114.  26518 

* 

-179. 82743 

* 

8 

7 

10 

H 

1.  11434 

* 

108.  80976 

58.  37181 

* 

1 

£ 

11 

H 

1.  1 1415 

* 

109.  01677 

* 

-57. £1596 

* 

1 

2 

1£ 

H 

1.  1 1 4£4 

* 

109.  02i98 

* 

122.  16535 

* 

£ 

3 

13 

H 

1.  1 14E6 

* 

108.  99666 

* 

- 1 2£.  20836 

* 

£ 

3 

14 

H 

1.  11419 

* 

109.  03092 

* 

57.  79603 

* 

3 

2 

15 

H 

1.  1 1 426 

* 

108. 98335 

* 

-57. 83985 

* 

3 

2 

16 

H 

1.11 4£4 

* 

108. 97829 

* 

57.  45893 

# 

4 

17 

H 

l.  1  14£9 

* 

108. 99140 

* 

. -58. 15307 

* 

4 

3 

18 

H 

1.  11423 

* 

109. 03421 

* 

57.  72829 

* 

5 

4 

19 

H 

1.  1  1426 

* 

108.  99207 

* 

-57.  90137 

* 

5 

4 

£0 

H 

1.11 423 

* 

109. 00945 

* 

57.  51033 

■* 

6 

5 

£1 

H 

1.  1 1428 

* 

109. 00320 

* 

-58. 12697 

* 

6 

5 

££ 

H 

1. 1 1419 

* 

109. 01732 

* 

57.  89713 

* 

7 

6 

£3 

H 

1.  11425 

* 

108. 99947 

* 

-57. 71635 

* 

7 

6 

£4 

H 

1.11 430 

* 

108. 96465 

* 

58.  1 1453 

* 

8 

7 

£5 

H 

1.  11410 

* 

109.  12599 

* 

-57. 59953 

* 

8 

7 

£6 

Tv 

10.  36744 

* 

. 00000 

.  00000 

1 

9 

NC 


1 

1 

1 

1 

£ 

£ 

3 

3 

4 

4 

5 

5 

6 
6 
8 


B-7 


(u  n  4-  m  ud  n  n 


INTERATOMIC  DISTANCES 


<^NOT  SHOWN 


EIGENVALUES 


-40.  61217 
-16. 77474 
-13. 41 163 
2. 99737 
4.  080 69 
5. 42985 


-36. 60451 
-15. 57225 
-12. 53984 

3.  £2518 

4.  37398 
5. 51490 


-36. 60281 
-15. 57103 
-12. 53792 
3. 48490 
4.  37472 
5. 51552 


-27. 76669 
-14. 32777 
-12. 27044 
3. 48524 
4.  65167 
5. 62519 


-27. 76660 
-14. 32688 
-12. 27004 
3. 79054 
4.  65202 
5. 62539 


-22. 10093 
-14. 25581 
-12. 08943 
3. 84729 
4. 69153 
5. 63487 


-21. 92546 
-14. 25422 
-1 1. 88432 
3. 84814 
4.  69217 
5. 63557 


-21.  92283 
-13. 41 169 
-11. 55576 
4. 08030 
4. 96231 
5. 64057 


NET  ATOMIC  CHARGES  AND  DIPOLE  CONTRIBUTIONS 


ATOM 

NO.  TYPE 

CHARGE 

ATOM  EL 

1 

C 

-.  0087 

4.  0087 

2 

C 

-. 0087 

4.  0087 

vJ 

C 

-.  0086 

4.  0086 

4 

C 

-.  0086 

4. 0086 

5 

C 

-.  0085 

4.  0085 

6 

c 

-.  0086 

4.  0086 

7 

c 

-.  0085 

4.  0085 

a 

c 

-. 0085 

4. 0085 

9 

H 

.  0041 

.9959 

10 

H 

.  0044 

.  9956 

1 1 

H 

.  0043 

.9957 

12 

H 

.  0042 

.  9958 

13 

H 

.  0043 

.9957 

1  4 

H 

.  0042 

.  9958 

15 

H 

.  0043 

.9957 

16 

H 

.  0043 

.  9957 

17 

H 

.  0043 

.  9957 

ia 

H 

.  0042 

.  9958 

19 

H 

.  0043 

.  9957 

20 

H 

.  0043 

.  9957 

21 

H 

.  0043 

.  9957 

C  C 

H 

.  0043 

.  9957 

23 

H 

.  0042 

.9958 

24 

H 

.  0045 

.  9955 

DIPOLE 

X 

Y 

Z 

TOTAL 

POINT- 

CHG. 

.  010 

.  006 

.  003 

.012  I — 

HYBRID 

.  001 

.  001 

002 

.  002  J 

SUM 

.  01  1 

.  007 

.  001 

.  013  «r - 

Spurious.  Due  to 
slight  charge  shifts. 
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CARTESIAN  COORDINATES 


NO. 

ATOM 

X 

Y 

Z 

1 

C 

.  0000 

.  0000 

.  0000 

2 

C 

1. 5435 

.  0000 

.  0000 

3 

C 

2.  1764 

1. 4078 

.  0000 

4 

c 

3.  7197 

1. 4077 

-.  0002 

5 

c 

4. 3528 

2. 8151 

-. 0089 

6 

c 

5. 8961 

2. 8149 

-.  0074 

7 

c 

6. 5294 

4. 2223 

-. 0238 

8 

c 

8.0729 

4. 2214 

-.0246 

9 

H 

3593 • 

.  5532 

-. 8982 

10 

H 

-.3630 

.  5704 

.  8856 

11 

H 

1. 9060 

-. 5612 

-. 8917 

12 

H 

1 . 9069 

-. 5611 

.8914 

13 

H 

1.8134 

1 . 9693 

.  8913 

14 

H 

1. 8135 

1.9683 

-. 8920 

15 

H 

4.0819 

.  8410 

-. 8886 

16 

H 

4.0824 

.  8517 

.8947 

17 

H 

3.9891 

3. 3827 

.  8783 

18 

H 

3.9908 

3. 3701 

-.  9048 

IS 

H 

6.  2597 

2. 2432 

-. 8920 

20 

H 

6. 2580 

2. 2644 

.8912 

21 

H 

6.  1670 

4.  7938 

.  8613 

22 

H 

6. 1659 

4.7732 

-.9216 

23 

H 

8. 4343 

3. 6544 

-. 9132 

24 

H 

8. 4381 

3. 6674 

.  8704 

ATOMIC  ORBITAL 

ELECTRON 

POPULATIONS 

1 . 22066 

.92761 

93138 

. 92900 

1. 22075 

.  92761 

.93130 

.92904 

1 . 22070 

.  92752 

93131 

. 92905 

1 . 22067 

.  92753 

. 93131 

. 92907 

1 . 22066 

. 92747 

93131 

. 92908 

1. 22069 

. 92748 

.93135 

. 929CG 

1.22069 

. 92738 

93133 

.92910 

1. 22063 

. 92743 

. 93135 

. 92904 

.  99586 

. 99559 

99572 

. 99578 

. 99567 

. 99578 

. 99572 

.  99571 

.  99568 

. 99575 

99572 

. 99571 

. 99573 

. 99568 

. 99582 

.  99553 
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APPENDIX  C 


USING  GIP  TO  MAKE  BIG  SYSTEMS 

The  GIP  function  of  the  DRAW  program  can  be  used  to  assemble 
two  or  more  systems,  A,  B,  C, . . . ,  to  form  the  single  system, 
A+B+C+ . . .  We  start  with  data  files  for  each  system  (let's 
say  there  are  three):  A. DAT,  B.DAT,  and  C.DAT. 

1.  We  will  select  system  A  to  be  at  the  origin  of  the 
cartesian  coordinate  description  of  the  A+B+C  system.  Bring 
A. DAT  into  DRAW  and  issue  the  GIP  command  as  follows:  GIP 
RETURN.  The  response  “OFFSET  IS  0"  will  be  printed  to 
indicate  that  the  coordinates  and  connectivity  will  be 
retained  as  in  the  original  A . DAT  file,  but  expressed  in 
cartesian  coordinates,  in  the  newly  created  A. GIP  file.  QUIT 
DRAW  without  performing  OUTPUT.  The  new  file  A. GIP  will  be 
in  your  directory. 

2.  Bring  B.DAT  into  DRAW  and  use  MOVE  to  place  it  (in 
cartesian  coordinates)  where  you  want  it  relative  to  the  A 
system  as  follows: 

In  DRAW> ,  issue  M0VE>. 

In  response,  type  CART  to  indicate  that  you  will 
designate  a  new  position  for  system  B  in  cartesian 
coordinates . 

You  will  be  prompted  for  information  on  where  you  want  B 
moved  to.  Note  that  B  can  also  be  rotated. 

After  the  desired  position  and  orientation  of  B  in  the 
cartesian  coordinate  system  (centered  on  A)  are  achieved, 
enter  GIP  #,  where  #  is  the  number  of  the  last  atom  in  the 

A.  DAT  file.  Hitting  RETURN  will  give  a  message  that  this 

B. GIP  file  is  “offset  by  Again  QUIT  DRAW  without 
performing  OUTPUT. 

3.  Bring  C.DAT  into  DRAW  and  use  MOVE  to  position  it 
relative  to  A  as  you  did  for  B.  Then  enter  GIP  #,  where  #  is 
the  sum  of  the  atoms  in  the  A. DAT  and  B.DAT  files.  This 
creates  the  C.GIP  file  "offset  by  #"  atoms.  QUIT  DRAW 
without  doing  OUTPUT. 

4.  Now  call  A. GIP  into  your  editor  (On  the  VAX,  EDT  A. GIP). 
At  the  end  of  the  A. GIP  file,  add  B.GIP.  (On  the  VAX,  place 
the  cursor  at  the  end  of  file,  hit  GOLD  7  and  on  prompt, 
enter  INCLUDE  B.GIP  and  hit  7  again.)  Then  delete  the  lines 
between  the  end  of  the  geometry  definition  in  the  A. GIP  file 
and  the  start  of  the  geometry  definition  of  the  B.GIP  file. 

5.  Add  the  C.GIP  file  to  the  now  combined  A  and  B  files  as 
described  above.  What  you  have  now  is  an  ABC. GIP  file  with 
cartesian  coordinates  in  the  geometry  definition  but  a 
connectivity  list  for  internal  coordinates.  Exit  the  editor 


renaming  the  combined  file  ABC.GIP  or  whatever  is 
appropriate . 

6.  Bring  ABC.GIP  into  DRAW  and  check  that  the  arrangement  of 
A,B,  and  C  is  as  desired.  The  EDIT  function  of  DRAW  may  be 
useful  in  making  needed  adjustments.  Then  issue  the  OUTPUT 
command  to  change  ABC.GIP  into  an  ABC. DAT  file.  The  ABC . DAT 
file  will  appear  in  your  directory  with  internal  coodinates 
and  the  proper  connectivity  list.  This  ABC. DAT  file  can  now 
be  operated  on  by  MOPAC . 


02 


